{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CSDN csdn_inside: [GMM 高斯混合模型 聚类 Python实现](https://blog.csdn.net/csdn_inside/article/details/85267341)\n",
    "\n",
    "[西瓜书数据集4.0](https://blog.csdn.net/xxliu_csdn/article/details/86483816)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "原理也不多说了。。大致思路就是把数据建立成k个高斯分布，EM迭代N次。最后看每个点在哪个高斯分布的概率最高，就分到那个分布。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里的computeGamma函数，用来算第i个簇的后验概率γji，用了2层循环，效率不高，本来想向量化的，搞了半天没搞出来。。干脆就先循环吧。。包括后面的fit方法也是一样。。用了很多循环。\n",
    "\n",
    "X的shape是(n_samples,n_features)\n",
    "\n",
    "mu的shape是(n_clusters,n_features)\n",
    "\n",
    "sigma的shape是(n_clusters,n_features,n_features)\n",
    "\n",
    "alpha的shape是(n_clusters,)\n",
    "\n",
    "输出的gamma是(n_samples,n_clusters)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "类方法fit里面，所有alpha初始为1/n_clusters，所有sigma协方差初始化为0.1，mu随机在样本里面选n_clusters个。\n",
    "\n",
    "这里为了复现西瓜书的效果，所以指定mu为[[.403,.237],[.714,.346],[.532,.472]]。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-01-29T04:34:53.682497Z",
     "start_time": "2020-01-29T04:34:52.252298Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    " \n",
    "def multiGaussian(x,mu,sigma):\n",
    "    return 1/((2*np.pi)*pow(np.linalg.det(sigma),0.5))*np.exp(-0.5*(x-mu).dot(np.linalg.pinv(sigma)).dot((x-mu).T))\n",
    " \n",
    "def computeGamma(X,mu,sigma,alpha,multiGaussian):\n",
    "    n_samples=X.shape[0]\n",
    "    n_clusters=len(alpha)\n",
    "    gamma=np.zeros((n_samples,n_clusters))\n",
    "    p=np.zeros(n_clusters)\n",
    "    g=np.zeros(n_clusters)\n",
    "    for i in range(n_samples):\n",
    "        for j in range(n_clusters):\n",
    "            p[j]=multiGaussian(X[i],mu[j],sigma[j])\n",
    "            g[j]=alpha[j]*p[j]\n",
    "        for k in range(n_clusters):\n",
    "            gamma[i,k]=g[k]/np.sum(g)\n",
    "    return gamma\n",
    " \n",
    "class MyGMM():\n",
    "    def __init__(self,n_clusters,ITER=50):\n",
    "        self.n_clusters=n_clusters\n",
    "        self.ITER=ITER\n",
    "        self.mu=0\n",
    "        self.sigma=0\n",
    "        self.alpha=0\n",
    "      \n",
    "    def fit(self,data):\n",
    "        n_samples=data.shape[0]\n",
    "        n_features=data.shape[1]\n",
    "        '''\n",
    "        mu=data[np.random.choice(range(n_samples),self.n_clusters)]\n",
    "        '''\n",
    "        alpha=np.ones(self.n_clusters)/self.n_clusters\n",
    "        \n",
    "        mu=np.array([[.403,.237],[.714,.346],[.532,.472]])\n",
    "        \n",
    "        sigma=np.full((self.n_clusters,n_features,n_features),np.diag(np.full(n_features,0.1)))\n",
    "        \n",
    "        for i in range(self.ITER):\n",
    "            gamma=computeGamma(data,mu,sigma,alpha,multiGaussian)\n",
    "            alpha=np.sum(gamma,axis=0)/n_samples\n",
    "            for i in range(self.n_clusters):\n",
    "                mu[i]=np.sum(data*gamma[:,i].reshape((n_samples,1)),axis=0)/np.sum(gamma,axis=0)[i]\n",
    "                sigma[i]=0\n",
    "                for j in range(n_samples):\n",
    "                    sigma[i]+=(data[j].reshape((1,n_features))-mu[i]).T.dot((data[j]-mu[i]).reshape((1,n_features)))*gamma[j,i]\n",
    "                sigma[i]=sigma[i]/np.sum(gamma,axis=0)[i]\n",
    "        self.mu=mu\n",
    "        self.sigma=sigma\n",
    "        self.alpha=alpha\n",
    "        \n",
    "    def predict(self,data):\n",
    "        pred=computeGamma(data,self.mu,self.sigma,self.alpha,multiGaussian)\n",
    "        cluster_results=np.argmax(pred,axis=1)\n",
    "        return cluster_results\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-01-29T10:01:58.049556Z",
     "start_time": "2020-01-29T10:01:58.039584Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "data=pd.read_csv('data_watermelon/watermelon_4.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-01-29T10:01:59.256367Z",
     "start_time": "2020-01-29T10:01:59.253336Z"
    }
   },
   "outputs": [],
   "source": [
    "# X = np.array(data.values, dtype=float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-01-29T10:49:21.243358Z",
     "start_time": "2020-01-29T10:49:19.975591Z"
    }
   },
   "outputs": [],
   "source": [
    "X = data.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-01-29T10:49:22.481037Z",
     "start_time": "2020-01-29T10:49:21.956410Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x2a6c4e45908>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VPW9//HXZ2aSycISIEFkUXYEwTWixX2pRa3grlS9Wut2LVeuXa71ahWt91ptq9KW21/dbV1wV1RcKu6KSlQEBFF2IlvYlyQzmZnv74+AhmQgA2TmzJy8n48HDzNnTua8j0PenHzPme8x5xwiIuIvAa8DiIhIy1O5i4j4kMpdRMSHVO4iIj6kchcR8SGVu4iID6ncRUR8SOUuIuJDKncRER8KebXh0tJS17NnT682LyKSkz799NNVzrmy5tbzrNx79uxJRUWFV5sXEclJZrYolfU0LCMi4kMqdxERH1K5i4j4kMpdRMSHVO4iIj6kchcR8aGUyt3MhpvZHDOba2a/SfL8xWZWZWbTtvy5tOWjiohIqpq9zt3MgsB44IdAJTDVzCY652Y1WvUJ59zoNGQUEZGdlMqR+1BgrnNuvnMuCkwARqY3loiI7I5Uyr0bsKTB48otyxo708ymm9nTZtYj2QuZ2eVmVmFmFVVVVbsQV0REUpFKuVuSZa7R4xeBns65/YA3gIeTvZBz7h7nXLlzrrysrNmpEUREssLXq1fxf1M/5t7PplK5Yb3XcVKSytwylUDDI/HuwNKGKzjnVjd4eC9w++5HExHx3h8/fJ8Hpn1KXTxO0ALcOeVDxh59HOcOHuJ1tB1K5ch9KtDPzHqZWT5wHjCx4QpmtmeDhyOA2S0XUUTEGzNXruCBaZ9SG4sRd45oIk4kHmPsO5Opqt7sdbwdarbcnXMxYDTwGvWl/aRz7kszu8XMRmxZ7Woz+9LMvgCuBi5OV2ARkUx5+es5RGOxJssDZrw5f54HiVKX0pS/zrlJwKRGy25s8PV1wHUtG01ExFtmRv1px21PMxoGlux0ZPbQJ1RFRLbjx/0HkB8KNlkedwmO79XHg0SpU7mLiGzHoLLOXHHQIYSDIfICAcLBIOFgkFuPPYHSoiKv4+2QZ3diEvED5xxEP4bYXAj1gfzDtvwqL34x5rBhnDpgHyYvmEcoEGR4n37s2bat17GapXIX2UUusQG35gKILwYXBwtCsAd0fAQLtPc6nrSg3h060rtDR69j7BQNy4jsIrfxfyA2D1w1EKn/b2w+bsOtXkcTUbmL7LKaSUBdo4V1UPuKF2lEtqFyF9llie0sj9ePxYt4SGPukhW2lmFOnYwMHwmRt9m25AOQf2Ru7UcOWVVdzeQF8zDguF59sv6KFS+p3MVTru4r3IaxUPc5EMYVnoG1uxazQq+jNcva3YRbPR0S1UA1UAiBIqz9TV5H86Unv5zBTW9PJmj1Aw43vT2ZW489gTMHDfY4WXZSuYtnXHw5bs0ocFvn6KiFmmdw8SVYx/s9zZYKC+4Jpf/C1bwEsdkQ2gcr/DEWaON1NN/5duMGbnp7MpF4HIh/t/yGt95gWI+9c+LSxExTuYtnXPU/wTU+IRmB6FRcbD4W6u1Jrp1hgWKs+FyvY/jeK9983WSe8a1enfcNPz3goIzmyQU6oSreqZsNRJsutxDE5mc8jmSvukScRJKT1AnnqIvHk3yHqNzFO3mDgfymy10dhPpmPI5kr+N79SEUaFpXAQtwQu/snuPFKyp38YwVXQAWZtubfYUhPAwL9fQolWSj/p1KueSAgygIhQgAAYyCUIjLDy7fqU+ObopGmbN6FRsikfSFzRIacxfPWLAzdHqi/hOd0U/ACqHwbKztL7yOJlnoV8OO5Ed9+/Py119hGD/uP4B9O++R0vcmnOP299/lH9OnEQoEiCXinLvvEH571LEEk/xG4Acqd/GUhfpiHR/yOobkiCGd92BIioXe0L2fVfDIjGlE4jEiW4bon5o1k5KCQv7zsGEtnDI7+POfLBGRBu77rIKaRndUqonFeHDaZx4lSj+Vu/iCq5tFYs0VJFYeSWL1BbjIx15HkiyyPlKbdPnGaCTpVTh+oHKXnOeiX+BWnwfRtyGxAuo+wa29jETNv7yOJlliYGlZ0uX9OnYi4NOpIlTukvPcxt8DtWx7n8ta2HirJvASAG48+lgKQqHvrssyoCAUYuzRx3kZK610QlVyX2x28uWJlfVzrFtxZvNI1jl4z248c/Yo/vzJFGZXVdG/Uymjhx7Gfnt08Tpa2qjcJfcFOkG8uulyC4MVZD6PZKWBZZ352ykjvY6RMRqWkdxXfCXQeBbJAig8H7Omd64XaQ105C45zwrPwiXWwub/q1/g4lB0Dtb2Gm+DiXhI5S45z8ywNpfjii+G+HIIlGIB3cRBWjeVu/iGWT6E9vI6hkhW0Ji7iIgPqdzF/xpf665r36UVULmLv40dC9dc832hO1f/eOxYL1OJpJ3KXfzLOVi3DsaN+77gr7mm/vG6dTqCF1/TCVXxLzO46676r8eNq/8DMGZM/XKfzikiAmBezb1RXl7uKioqPNm2tDLOQcMbMiQSKnbJWWb2qXOuvLn1UhqWMbPhZjbHzOaa2W92sN5ZZubMrNkNi2TE1qGYhhqOwYv4VLPlbvWf3x4PnAQMAkaZ2aAk67UFrgY0kbZkh4Zj7GPG1B+xjxmz7Ri8iE+lMuY+FJjrnJsPYGYTgJHArEbr/Q64A/hViyYU2VVmUFKy7Rj71jH4khINzYivpVLu3YAlDR5XAoc2XMHMDgR6OOdeMjOVu2SPsWPrj9C3FvnWglexi8+lMuae7Kfgu99nzSwA3AX8stkXMrvczCrMrKKqqir1lCK7o3GRq9ilFUil3CuBHg0edweWNnjcFhgMvG1mC4HDgInJTqo65+5xzpU758rLypLf9kpERHZfKuU+FehnZr3MLB84D5i49Unn3HrnXKlzrqdzrifwETDCOafrHEVEPNJsuTvnYsBo4DVgNvCkc+5LM7vFzEakO6CIiOy8lD6h6pybBExqtOzG7ax7zO7HEhGR3aG5ZUREfEhzy0hOc7G5EJkCgfYQPh4LFHsdSQSAGStX8NmybykrasPxvXoTDmW2blXukpOcc7gNN0HN84ADCwJjocP9WP6BHqeT1iyWSHDVyxP5YMkiEs4RCgQIh0JMOPNc+nbslLEcGpaR3BSZDLUvALVABFw1uE24tVfiXNzrdNKKPT5zOh8sWURNLEYkHmdzXR1ra2q46uWJzX9zC1K557hopI77rnuUszpfwqntLuSWc/7EysX+/4CYq34KXE2SZ6JQ93n6tutiuNo3SGy8E1f9BC6xKW3bktw0YeZ0amKxbZY5oHLjBhavX5exHBqWyXE3nXY709+ZRbS2DoAPnv2YGe/M4oGvxtG2QxuP06VTdDvLDVxsO8/tHpfYhFszCuJLwFXjKISNf4COj2N5/dKyTck9dYnkvzkaRiyRyFgOHbnnsAUzFjHjvdnfFTtAIuGo2VTLqw+86WGy9LPC04DC5E/mH5SWbbrNf4PYgvohIABqwG3ErW925g1pRU4bMIhwsOlxc4eCAnqVdMhYDpV7DlswYzGBYNO3MFIT5auPv/EgUQYVnAL5h4AVbVmQDxRg7f9E/Qep06DmRZr+xuAgNg+XWJOebUrO+ekBB9G/UyeK8vIACAdDFOflMe6kU7AMzmukYZkc1rVvF1yi6Zzk+QV59By8lweJMscsBB3ugegUXOQ9CHTACkdgwT3TuNUdHQtpMjKpV5iXxzPn/IS3Fszn428r6dq2LSMHDKRTUVHz39yCVO45bMAhfdlrYDfmT19MLPr9OHMoP8Qpl5/gYbLMMAtA+HAsfHhmNlh4Omy+D4g0WBiAvIFYIHO/bkv2CwUC/LBPX37Yp69nGTQsk8PMjNtfv5EjzjiUUH6IQDDAgKF9ufOdW+jYRWXT0qzNFZA3aMtQUAisGAIdsfZ3eh1NpAndINsnYnUxEvEE+QVpGm8WoP7DU0Q/hthMCHSFghPSN8YvkkSqN8jWsIxPhPJCkOd1Cv8zMwgfVv9HJItpWEZExId05C4ivvLCnNnc/dGHLN+0iX6dOnHd4Ufxgx7+vnosGR25i4hvPDp9Gv89+XUWrV9HJB5j5soV/OzF5/ioconX0TJO5S4ivpBwjj999EGTeV1qYzHu+OBdj1J5R+UuIr6wMRJhczT5nENz17a+TxCr3EXEF9rk5yed0wWge7t2GU7jPZW7iPhCMBDgivJDKGx0x6OCUIhfHJahTzFnEV0tIyK+cVX5oQQtwP+r+ITNdVFKi4q47oijOaG3d9MAeEXlLiK+YWZcWT6Uyw8+hGg8RjgYyuhMjNlE5S4ivhMwoyDUuj+yrTF3EREfUrmLiPiQyl2klVhVXc17ixYyd81qr6NIBmjMXcTnnHPc+t7bPDbjC/KDQeoSCQaWlnH/iNMpKdjOfWgl5+nIXcTnnp79JRNmTicSj7MxGqU2Vj/nyjWvveJ1NEkjHbmL77joJ7hNf4f4Msg/FGtzBRbs4nUszzzw+adN5lupSySYUrmYdbU1Onr3KZW7+Eqi+jnYcBNQW7+gZiGu9kUofQELdvM0m1fWR2qTLg+YsSkaVbn7lIZlxDecq4ON/8N3xQ5ADNxm3KbxXsXy3DF79yKU5IM8bfPDdG3b+uZcaS1U7uIf8UogluwJiHyY6TRZY8yhw2hfUEg4GAQgaEZhKMRtx59IoJV+erM10LCM+EegBFyycgeCZZnNkkX2aNOG1y64iH9On8aUJYvZq30Jlxx4MPuUtt7/J61BSuVuZsOBcUAQuM859/tGz18J/ByIA5uAy51zs1o4q8gOWaADLnwURN4FGszrbYVY8eWe5coGHQuLGHPoMMYcOszrKJIhzQ7LmFkQGA+cBAwCRpnZoEarPeacG+KcOwC4A7izxZOKpMDa3wH5PwDCYG2AQij+D6zgh15HE8moVI7chwJznXPzAcxsAjAS+O7I3Dm3ocH6xYBryZAiqbJAG6zjvbj4CkisglBvzHQ1iLQ+qZR7N6Dh3WUrgUMbr2RmPwd+AeQDx7VIOpFdZME9ILiH1zFEPJPK1TLJTqc3OTJ3zo13zvUBrgVuSPpCZpebWYWZVVRVVe1cUhERSVkq5V4J9GjwuDuwdAfrTwBOS/aEc+4e51y5c668rExn6kVE0iWVcp8K9DOzXmaWD5wHTGy4gpn1a/DwFOCblosoIiI7q9kxd+dczMxGA69RfynkA865L83sFqDCOTcRGG1mJwB1wFrgonSGFhGRHUvpOnfn3CRgUqNlNzb4ekwL5xIRkd2g6QdERHxI5S4i4kMqdxERH1K5i4j4kC/KPZFIEI3UeR1DRCRr5HS5RyN1jL/6fka0vZAfF5/PpYOv4Yt3vvQ6loiI53K63P/w0/FMuv9NIjVRXMKxaFYl159yGwtmLvY6mogvxBMJHp/xBSMnPMLJj/2Dez6dSm1MvyXngpwt99XL1vLB858QrYlus7yuNsoTtz/vUSoRf/mPV17i1vfeZsbKFXy1qoq7P/6QUc88SSyR8DqaNCNny335gpXkh/OaLE8kHAtnLknyHSKyM2ZVreSdRQuoiX1/d6vaWIxv1qzmrQXzPUwmqcjZcu/ef0/qkpxEDYYCDBjad7dee/2qDfxl9H2c2/UyLuzzc56443liddu5fZuIT1Us/ZaEa3prhuq6OqZUpmnos/H2kmxfUpOz5d6+tB0nXnws4aL8bZbnF+Rz7n+N3OXXrdlcy1WHXMuke99gzfJ1LF+wkn/e/BS/O1c3l5LWpbSomFAg2GR5OBikS5u2Lb/BsWPhmmu+L3Tn6h+PHdvy22oFcrbcAUb/5RIuvPFsOu7ZgfyCfA44bjB3v38rXft02eXXnPzIe2xYtZFYXfy7ZZGaKBWvfaETtdKqHN+rN/nBYJMbOgQDAU4f2PhOm7vJOVi3DsaN+77gr7mm/vG6dTqC3wUpTRyWrYLBIOf+12mc+19Jp4/fJTPfn03t5kiT5YGA8c2n8+k1eK8W25ZINguHQkw481yuePl5VmzahJnRJj+fPw//MWVFxS27MTO46676r8eNq/8DMGZM/XJLds8g2ZGcLvd06NZvT/LCeU3G882MPfbWDUakdenXqROTL7yEBevWEksk6NuxE4F0Fe3Wgt9a7KBi3w05PSyTDiddejyhvG3HGYOhAJ26dmDIUQM9SiXiHTOjd4eO9O9Umr5ih++HYhpqOAYvO0Xl3khp147c8caNdO/flbxwiFB+iCFHDeKPb91MIJC5/10bVm9kTsU8NqzZmLFtinim4Rj7mDGQSNT/t+EYvOwUDcsksc/Qfjz41TjWrlhHXjiPNiUtPL64A/FYnHFX3cvkR94llB8iFo1x4kXHMPqvPyMYbHrlgogvmEFJybZj7FvH4EtKNDSzC8x59C9ieXm5q6io8GTb2ezB307gmbteJFL9/Sdvw0X5nHftaVzw27M9yTR32gLeffojAgHj6HOG6aSypI9z2xZ548eCmX3qnCtvdj2Ve3Y5rcNFbF5f3WR5u05teKbqwYznefCGx3nm7peoq60DM/LyQ5z/27MY9ZvTM55FRFIvd425ZxHnHNUbapI+t2ld08JPtwUzF/PMXS8RqY6SSDgS8QSRmiiP3PIUy+avyHgekUzbHI3y3OxZ3PdZBV+sWO51nJ2iMfcsYmb0OWBv5n6+sMlz/Q7unfE8Hzz/CXXRptMuOAdTJlZwxn+ekvFMIpkyc+UKzn/2KRIuQTQeJxQIcOTePRl/0qkEM3hxxa7K/oStzM///DPCRWECgfpxxkDACBeF+fm4SzKeJRQKfpejIQsYgZD+6oh/Oee44qUX2BiNsLmujrpEgppYjPcWLeSZ2blxzwj9hGaZwYfvw18++l+OOmcYew/qzjHnHcFfP76NgYf2y3iWo87+AYFQkit0nOOIMw7NeB6RTJmzehXrI7VNltfEYjzx5QwPEu08DctkoV6D9+L6x/7T6xh07dOFK/5wIX//1T8wMzDDJRJc/bfLKO3a0et4ImmTcK7JnDpbxXPkmnuVu+zQiKuGM2zkIUx58VMCAWPYyEPosEeJ17FE0mqf0jIK8/LYXLftNCSFoRBnDdzXo1Q7R+UuzSrt1olTrzzR6xgiGRMw468nncolE58l4Ry1sRhFeXns17kL5+w7xOt4KVG5+5U+DJLVqipXM+2tmRS3L6L8RwckvauYeGtot+68e/GlTJzzFVXVmzmsWw8O32vv9M6v04JU7n40dmz9HNhbP8a9dd6OkhLd+CALPHTTBJ78w0SCoSABM4J5QX7/2g30P7iP19GkkY6FRVx8wEFex9glulrGb7L4pgeRmgjP/WUSvzzmJm46/Q4+e2O6Z1m88vmbM3jmzvpP/NZuqqV6Yw0b12zi+lNuIx6PN/8CIinSkbvfZOlND6K1Ua4edj3ffrPsu3lzPv3XdEZddzrnX3+mJ5m88PLf/5X0ZjDRmigz3/+K/Y/OjZN1kv105O5HDQt+K49vevDGI++xdO7ybSZEi1RHePTWZ1hXtd6zXJlWvanptdMAWH3Bi7QUlbsfZeFND6ZMnJr0iDUvHGLWh197kMgbx553OAXF4SbL47EEg4/Yx4NE4lcqd7/J0pselOzRPulUBs452nZs40Eibxx73uEMOKQvBcUFQP1dvsKF+Yz522UUtin0OJ34SUpj7mY2HBgHBIH7nHO/b/T8L4BLgRhQBVzinFvUwlklFVl604NTrzyRtx57n0iDoQczaFPShn0PH+BJJi+E8kLc/vpvmfJiBR9OnEq7Tm056ZLj2HtQD6+jtagNkQifL1tK23CYA7rsmTOXD/pJs/O5m1kQ+Br4IVAJTAVGOedmNVjnWOBj51y1mf07cIxz7twdva7mc0+zLLzO/ZUHJjP+6gcJhgK4hKOkczv+d9L1dO/f1dNc0rIe/uIzfv/+e+QHAySco31BAQ+PPJM+HTt5Hc0XWuxmHWb2A2Csc+5HWx5fB+Ccu2076x8I/NU5d/iOXlfl3jrVbK5lzidzKWpXSL+DetfPWSO+8dmypVz43FPUxL6fKtqALm3a8t5PL9MRfAtoyZt1dAOWNHhcuWXZ9vwMeGU7oS43swozq6iqqkph0+I3hcUFHHDsYPof3EfF7kOPTJ9GbWzbewA4YEOkls+WLfUmVCuVSrkn+wlMerhvZhcA5cAfkj3vnLvHOVfunCsvKytLPaXsUF20jmULVlCzeTuX2YlkyJramqTlYGZsiDS9WkrSJ5UTqpVAw7M93YEm/wSb2QnA9cDRzjm9ixny9J0v8s+bnyKRSJBIOE6+9Hiu/NNFBJPNwy6SZj/q3Zep31ZuMywDUBePU95V51YyKZUj96lAPzPrZWb5wHnAxIYrbBln/zswwjm3suVjSjL/+uc7PHTjE1RvrKF2c4RoTZRX7p/M/f/9qNfRpJU6Y+C+9O7QkcJQ/XGjUT9N7q+GHUm7cIG34VqZZk+oApjZycDd1F8K+YBz7n/M7Bagwjk30czeAIYAy7Z8y2Ln3IgdvaZOqO6+SwaOYcmcpuOY4aIwz699iFCeZpeQzIvEYjw7+0tenfcNHQoKOX+//Tmka3evY/lGqidUU/rpd85NAiY1WnZjg69P2OmEsttWL1ubdHm8Lk7t5ghtSlTuknnhUIhRQ/Zn1JD9vY7SqukTqjms38G9ky5vV9qW4vZFGU4jItlE5Z7DLrv9QsJF4W0+mxQuyuff77pIlxmKtHIq9xw2oLwPd7//O4aecjCdunZg8BH7cMvz13LMOTv8/JiItAIpnVBNB51QFRHZeS35CVUREckxKncRER9SuYuI+JAuhBaRlCzftJE/fvg+by9cQFF+HhfudyCXHHAQwYCOEbORyl1EmrW+tpYREx5hbU0NcedYU1vDXR99wKyVK7hr+Clex5Mk9E+uiDTr8ZnT2RSJEm9wdV1tLMar875h8fp1HiaT7VG5i0izpi6tpDYea7I8Lxhk9irdmyEbqdxFpFm9O3QkL8nYejzh6Na2nQeJpDkqd/GNRCLBC+Nf4aJ+ozmz7BJuPe8uli1Y4XUsX/i3/Q4kFNj2HgF5gQB9O3Zk37LOHqWSHVG5i2+MH/MA9177KEvnrWDD6o289/QUriq/druzZ0rqerRvz8OnnUnPkg7kBYLkBQIc3bMXD592puYxylK6WkZ8Ye3K9bxy35vUReq+W5ZIOCKbIzx798tcdvsFHqbzh/Ku3Zh84U9ZU1NDQShEcX6+15FkB1Tu4gsLZiwmvyBvm3IHqIvGmPn+bI9S+Y+Z0alI00nnAg3LiC906VnWpNgBAsEAPQbo3p3S+qjcxRe69unCvofvQ154219G88J5nPmLUz1KJeIdlbv4xthnf82w04aSFw6RFw6xR88ybn7u1/QavJfX0UQyTmPu4htFbQu54fFrqK2OULu5lval7XQlh7RaKnfxnYKiMAVFYa9jiHhKwzIiIj6kchcR8SGVu4iID6ncRUR8SOUuIuJDKncRER9SuYuI+JDKXUTEh1TuIiI+pHIXEfEhTT8gIr713uKFPDztc9bW1jC8Tz9+MmT/VnOTkZSO3M1suJnNMbO5ZvabJM8fZWafmVnMzM5q+Zitl3OOzRuqidU1vfO8iGzf36Z+zJUvvcCbC+fz+fJl3PnRh5z2xKNU1zWd99+Pmi13MwsC44GTgEHAKDMb1Gi1xcDFwGMtHbA1++SVz/m3PqM5s/QSRpZcxF9G30c0yQ0pRGRb62pr+PMnU6iJfX9QFInHWLpxA0/PmulhssxJ5ch9KDDXOTffORcFJgAjG67gnFvonJsOJNKQsVWaM3Uut5z9R5YvXEk8FidaE+XVB9/izkv/5nU0kaz3+fJl5AWDTZbXxGK8MX+eB4kyL5Vy7wYsafC4cssySaPHb3uWaE10m2XRmijvPv0R66rWe5RKJDd0LCgk4VyT5QEzOhcXe5Ao81I5oZrsbgdN/6+l8kJmlwOXA+y1l+6OsyNL5iwlyd9N8sIhVlWuoaSsfeZDieyiFZs2cd3k1/lkaSWFoTx+esBBXHXIoWnb3n57dKGsqJglG9ZvU/L5wSAX7ndA2rabTVI5cq8EejR43B1Yuisbc87d45wrd86Vl5WV7cpLtBoDDulLINj07YlFY+zZZw8PEonsmlXV1Rz10L28vWgB1XV1rK6p5o9T3uei559O2zbNjH+cdha9SjpQGMqjTX4+haE8bjrqOPbvsmfatptNUjlynwr0M7NewLfAecBP0ppK+Ml/n8H7z35Mzaba75aFi8KMHD2c4nZFHiYT2Tk3vvUv6hJNT8e9t3gRi9atZe+SDmnZbo/27Xn9gouZs3oVGyIRhnTeg8K8vLRsKxs1e+TunIsBo4HXgNnAk865L83sFjMbAWBmh5hZJXA28Hcz+zKdoVuD7v27Mu6DWyk/cX8K2xbSpWdnLv/DhVx62/leRxPZKVMql2z3uRfmfJXWbZsZ+5SWMbRb91ZV7JDih5icc5OASY2W3djg66nUD9dIC+o1ZG9ue/UGr2OI7JZ2+WHWRyJJn+vatl2G07Qemn5ARNJq9NDDki4PBQKcsc/ADKdpPVTuIpJWZ+87hDMH7rvNsrxAgEdPP5tAQBWULuaSXW+XAeXl5a6iosKTbYtI5q2prubFb+bQubiYH/Xuq2LfRWb2qXOuvLn1NHGYiGREx6IiLtr/QK9jtBr6p1NExIdU7iIiPqRyFxHxIZW7iIgPqdxFRHzIF1fLrKtaz+sPv8OyecsZfMRAjjzrMPLDreujxiIiDeX8de5ffzqPXx93M7EtN7QoaFNAadcO/OWj22hT0jrmbRaR1iPV69xzfljm9xf+meqNNd/d2KJ2Uy3LF1bxyO/SN52oiEi2y+lyX7N8LcsXVDVZHovGeOepDz1IJCKSHXK63IOhINsbVsrL15i7iLReOV3u7Uvb0b+8T5M7FoUL8zn50uM9SiUi4r2cLneA6x8bQ2m3jhS2LSRcmE+4KMyQowZx5i9+7HU0kYxxzvHF8mW8Pu8blm3c6HUcyQI5fylk573K+Mfcv1Lx+hesXLyK/uV9GFDex+tYIhmzcvMmLnzuab7duIH8yRsKAAAFAklEQVQARjQR55xBg7n5mOMxS3Z/e2kNcr7coX7s/dCTD/I6hognRk96iflr1xBvcP7pmdmz2L/Lnk3mUZfWI+eHZURas5WbNzFj5fJtih2gJlbHQ9M+8yiVZAOVu0gO2xSNErDkP8Ybt3PfUmkdVO4iOaxnSQeKk1z2mxcI8qM+/TxIJNlC5S6SwwJm3HHCcApCIYJbTp4WhEKUFRdxZflQj9OJl3xxQlWkNTumZy9ePO8C/jF9GpUb1nN4j705e9Bg2obDXkcTD6ncRXygT8dO3HyMPrgn39OwjIiID6ncRUR8SOUuIuJDKncRER9SuYuI+JDKXUTEh1TuIiI+pHIXEfEhlbuIiA+p3EVEfEjlLiLiQ+YaTfKfsQ2bVQGLPNl4epQCq7wOkSZ+3jfw9/75ed/A3/u3vX3b2zlX1tw3e1bufmNmFc65cq9zpIOf9w38vX9+3jfw9/7t7r5pWEZExIdU7iIiPqRybzn3eB0gjfy8b+Dv/fPzvoG/92+39k1j7iIiPqQjdxERH1K57yQzG25mc8xsrpn9JsnzV5rZDDObZmbvm9kgL3Luiub2rcF6Z5mZM7OcuUohhfftYjOr2vK+TTOzS73IuatSee/M7Bwzm2VmX5rZY5nOuKtSeO/uavC+fW1m67zIuatS2L+9zOwtM/vczKab2ckpvbBzTn9S/AMEgXlAbyAf+AIY1Giddg2+HgG86nXultq3Leu1Bd4FPgLKvc7dgu/bxcBfvc6axv3rB3wOdNjyuLPXuVtq3xqt/x/AA17nbuH37h7g37d8PQhYmMpr68h95wwF5jrn5jvnosAEYGTDFZxzGxo8LAZy5aRGs/u2xe+AO4DaTIbbTanuW65KZf8uA8Y759YCOOdWZjjjrtrZ924U8HhGkrWMVPbPAe22fN0eWJrKC6vcd043YEmDx5Vblm3DzH5uZvOoL8GrM5RtdzW7b2Z2INDDOfdSJoO1gJTeN+DMLb/2Pm1mPTITrUWksn/9gf5m9oGZfWRmwzOWbvek+t5hZnsDvYA3M5CrpaSyf2OBC8ysEphE/W8nzVK57xxLsqzJkblzbrxzrg9wLXBD2lO1jB3um5kFgLuAX2YsUctJ5X17EejpnNsPeAN4OO2pWk4q+xeifmjmGOqPbu8zs5I052oJKf3MbXEe8LRzLp7GPC0tlf0bBTzknOsOnAz8c8vP4w6p3HdOJdDwiK47O/4VaQJwWloTtZzm9q0tMBh428wWAocBE3PkpGqz75tzbrVzLrLl4b3AwRnK1hJS+XtZCbzgnKtzzi0A5lBf9tluZ37mziO3hmQgtf37GfAkgHNuClBA/bwzO6Ry3zlTgX5m1svM8qn/yzSx4Qpm1vAH5hTgmwzm2x073Dfn3HrnXKlzrqdzrif1J1RHOOcqvIm7U1J53/Zs8HAEMDuD+XZXs/sHPA8cC2BmpdQP08zPaMpdk8q+YWYDgA7AlAzn212p7N9i4HgAMxtIfblXNffCoRYO6mvOuZiZjQZeo/4s9wPOuS/N7Bagwjk3ERhtZicAdcBa4CLvEqcuxX3LSSnu29VmNgKIAWuov3omJ6S4f68BJ5rZLCAO/No5t9q71KnZib+Xo4AJbsslJbkixf37JXCvmV1D/ZDNxanspz6hKiLiQxqWERHxIZW7iIgPqdxFRHxI5S4i4kMqdxERH1K5i4j4kMpdRMSHVO4iIj70/wE5t9O1Rm70GgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model1=MyGMM(3)\n",
    "model1.fit(X)\n",
    "result=model1.predict(X)\n",
    "plt.scatter(X[:,0],X[:,1],c=result)\n",
    "plt.scatter(model1.mu[:,0],model1.mu[:,1],marker='x',color='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "nbTranslate": {
   "displayLangs": [
    "*"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "en",
   "targetLang": "fr",
   "useGoogleTranslate": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
